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New stationary solutions of the (Michelson) Sivashinsky equation of premixed flames are obtained 
numerically in this paper. Some of these solutions, of the bicoalescent type recently described 
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conditions, the time evolution of the Sivashinsky equation in the presence of a moderate white 
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I. INTRODUCTION 

The Sivashinsky equation [l| (or Michelson Sivashinsky equation depending on the au- 
thors) is a well established non linear equation which provides a satisfactory description of 
the time evolution of premixed flames. Until working on the present paper, the author had 
a very simple idea of the situation regarding this equation. Pole solutions of the Sivashinsky 
equation were obtained in [2| and popularized in jj|, which reduces the time evolution of 
the equation to a dynamical system and the stationary solutions to finding zeroes of a non 
linear function of several variables. The paper |3| also shows that the poles have a tendency 
to coalesce, i.e. to align vertically in the complex plane. Stationary solutions were obtained 
in the form of coalescent solutions with a number of poles depending on the width of the 
domain. It was shown analytically in (j|p| that each solution, with a given number of poles 
is linearly stable in a given interval for the control parameter (either the domain width or 
more often the curvature term with a domain width fixed to 27r). Numerical simulations 
however, always performed with periodic boundary conditions, continue to show that the 
solutions are extremely sensitive to noise p| for sufficiently large domains. These results are 
consistent with a qualitative description of the stability of curved fiame fronts by Zeldovich 
et. al. y. 

For some reason, the author of the article began simulations of the Sivashinsky equation 
with Neumann boundary conditions, ie. zero slope of the fiame front at each end of the 
domain. Of course, Neumann boundary conditions are a more realistic description of a fiame 
in a tube than periodic boundary conditions. However, as solutions with Neumann boundary 
conditions on [0, vr] are simply symmetric solutions with periodic boundary conditions on 
[0, 2tt], the author was thinking that he should obtain basically a coalescent solution, but only 
between and vr, with all the poles coalescing at 0, leading to a cusp at this boundary. It was 
so obvious that actually simulations of the Sivashinsky equation with Neumann boundary 
conditions were only used originally as a test case for a new computer program. However 
stationary solutions were obtained, where poles did not all coalesce at the same position, 
but actually on the two boundaries. 

It turns out (although the author was absolutely unaware of this paper at the beginning 
of his work) that this type of stationary solutions, called bicoalescent solutions, were already 
discovered by Guidi and Marchetti p|. In Section HH we show the new bicoalescent solutions 



that we have obtained, which have a nice property with Neumann boundary conditions, they 
are stable. These solutions were not found in (s| because the curvature parameters studied 
were too large (or equivalently, the domain width was too small). In Section ITTH we show 
where the new solutions of Section Ull are found in theparameter space. We have thus to 

r 

study a larger domain of the parameter space than in p|, and discover also new stationary 
solutions of the interpolating type described by Guidi and Marchetti (see section IIIII for a 
definition of this type of solution). These interpolating solutions, unlike those of Section ITU 
are unstable. The number of stationary solutions obtained is so large that we have entitled 
Section Uni web of stationary solutions and will try to convince the reader that this is not an 
exageration. In Section HV| the evolution of the Sivashinsky equation with noise is studied. 
In the case of Neumann boundary conditions, as expected, the stable bicoalescent solutions 
play a dominating role in the dynamics. Finally, Section |3 contains a conclusion. 

II. STABLE BICOALESCENT SOLUTIONS 

The Sivashinsky equation can be written as 

<t>t + -(t>l = T^<Px. + I (0) (1) 

where (p (x) is the vertical position of the front. The Landau operator / (</>) corresponds 
to a multiplication by \k\ in Fourier space, where k is the wavevector, and physically to the 
destabilizing infiuence of gas expansion on the fiame front (known as the Darrieus-Landau 
instability), u is the only parameter of the equation and controls the stabilizing infiuence of 
curvature. The linear dispersion relation giving the growth rate a versus the wavevector is, 
including the two effects: 

a=\k\-uk^ (2) 

As usual with Sivashinsky-type equations, the only non linear term added to the equa- 
tion is |0^. In the fiame front case, this term is purely geometrical : the fiame propagates 
in the direction of its normal, a projection on the vertical (y) direction gives the factor 
cos (6*) = 1/a/1 + 0^, where 9 is the angle between the normal and the vertical direction, 
then a development valid for small slopes of the front leads to the term ^0^. The Sivashinsky 
equation will be solved numerically on [0, 27r] with periodic boundary conditions, or (more 



often in this paper) on [0, 27r] with only symmetric modes, which corresponds to homoge- 
neous Neumann boundary conditions on [0, tt] (zero slope on both ends of the domain). 
All dynamical calculations will be performed by Fourier pseudo-spectral methods (i.e. the 
non linear term is calculated in physical space and not by a convolution product in Fourier 
space). The method used is first order in time and semi-implicit (implicit on the linear terms 



of the equation, exr 
Pole solutions {[, 



icit on |0^). No particular treatment of aliasing errors has been used. 
) of the Sivashinsky equation are solutions of the form: 

,^2.|:{..(.„(i^)).,„(s,„(i^))) (3) 

where N is the number of poles Zn{t) in the complex plane. Actually the poles appear in 
complex conjugate pairs, and the asterisk in Equation [HI denotes the complex conjugate. In 
all the paper, only poles with a positive imaginary part will be shown, the number of poles 
will also mean number of poles with a positive imaginary part. The pole decomposition 
transforms the solution of the Sivashinsky equation into the solution of a dynamical system 
for the locations of the poles. In the case of stationary solutions, the locations of the poles 
are obtained by solving a non linear system: 

-^ E cot(^^—^)~tsgn[lm{zn)] = n=l,---,N (4) 

where Im (z„) denotes the imaginary part and sgn is the signum function. This non linear 
system will be solved by a Newton-Raphson method. 

Let us define here a process that will be called folding in the rest of the paper and which 
allows to create cellular solutions starting from curved flame fronts (i.e. fronts with only 
one cell in [0, 27r]). If a solution (pi (x) of the Sivashinsky equation exists with parameter 
1/z/i, then 02 (x) = —01 (mx) is a solution of the Sivashinsky equation with parameter 
l/z/2 = m (l/z/i),with m integer. 

Although we have searched for stationary solutions with periodic boundary conditions, it 
appears that all the solutions we have found on [0, 27r] are symmetric, and thus are stationary 
solutions with Neumann boundary conditions, i.e. zero slope, on [0,7r]. In most of the cases 
the stationary solutions obtained have poles at x = 0, in a few cases however, the solutions 
have no poles on the boundaries (i.e. only lead to symmetric solutions with no poles at the 
boundary) 



With periodic boundary condition, the well-known result is that in the window 2n — 1< 
l/z/< 2n + 1, n = 1,2, ■■ ■ there exists n different monocoalescent stationary solutions (all 
the poles have the same real part), with 1 to ra poles, and the solution with the maximum 
number of poles n is asymptotically stable. For a particular value of l/u, the number n(i/) 
such that 2n — l< l/z/< 2n+l will be called the optimal number of poles. All stable solutions 
found in this paper, for any value of 1/z/, even with Neumann boundary conditions, have 
the optimal number of poles n(z/). 

Using however the Sivashinsky equation (Eq. Q]) with Neumann boundary conditions, we 
obtain in each of the intervals [2n — l,2n + 1] of the parameter 1/z/, not only one asymp- 
totically stable solution, but several, of the form (/, n — I) with / = 0, 1,- ■ ■ , ra where / poles 
coalesce at x = and / — n coalesce at x = vr (The bicoalescent type of solutions have been 
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recently introduced in p|). These solutions will also be obtained from the non linear system 
of equations (Eq. ^ in Section IIIIl It must be remarked that all these solutions, except 
the monocoalescent one, are unstable for periodic boundary conditions, i.e. when antisym- 
metric perturbations are allowed on [0,27r]. We have just defined here the notation (ni,n2) 
that will be used in the paper for bicoalescent solutions with rii poles at zero, and n2 at vr. 
Monocoalescent solutions can be seen as a particular case of bicoalescent solutions and will 
be noted (n, 0). We will encounter also multicoalescent solutions, such as (ni,n2,n3, ■ ■ ■), 
which means that in the interval [0, 27r], the poles coalesce at different locations: rii poles 
coalesce at a position on the left of the interval, generally 0, n2 poles coalesce at a position 
with a higher value of x, then ns at a position with a value of x even higher, and so on. With 
this notation (1,1,1) represents a cellular solution with three cells obtained by the folding of 
the (1,0) solution. 

For the particular value 1/z/ = 10.5 (five poles) the different possible solutions are shown 
on [0, vr] in Figure [H On the left, we have a monocoalescent (5,0) solution with five poles at 
0. The middle solution of the figure is a (4,1) solution (4 poles at x = 0, 1 pole at x = vr). 
Finally the solution on the right is a (3,2) solution (3 poles at 0, 2 poles at vr). For an even 
value of the optimal number of poles (i.e. the value of n in the interval [2n — 1, 2n + 1]), 
the stable solutions will include a solution symmetric on [0, vr], for instance if ra = 6 we have 
the solutions (6,0) (5,1) (4,2) and the symmetric (3,3) solution. In Figure O, we show on 
the same figure the shape (x, 0(x)) (x is the horizontal direction) of the (3,2) solution for 
1/z/ = 10.5 with X G [0, vr] (lower part of the figure, below the horizontal segment) and the 



corresponding locations of poles in the complex plane (upper part of the figure, above the 
horizontal segment). The poles are indicated by circles, the segment is the real axis in pole 
space between and 27r. The important thing about this type of figure, which will be used 
in the rest of the paper, is that a pole very close to the real axis (i.e. very close to the 
horizontal segment in the upper part of the figure) leads to a cusp in the front shape (in the 
lower part of the figure), and that the x value of the pole in the complex plane is the same 
as the X value in physical space of the cusp that appears; in a diagram like Figure [2l the 
cusp and the corresponding pole are on the same vertical line . We will see later however 
examples of poles far away from the real axis with no cusp at the x value of the pole. This 
effect results from a competition between a new pole and the poles at zero which tend to 
prevent the appearance of a new cusp. It is described in a simple way in Appendix 1X1 

An illustration of the stability of the (3,2) stationary solution is given in Figure 01 The 
initial condition used in the Sivashinsky equation, with Neumann boundary conditions, is 
exactly the (3,2) solution for l/u = 10.5. In a simulation without noise, the amplitude 
(maximum minus minimum of (p{x)), would simply stay constant with time, as the (3,2) 
solution is stable. In order to complicate the convergence to the (3,2) solution, we apply a 
noise (additive gaussian white noise added to the Sivashinsky equation, amplitude a = 0.01, 
see section IIVI for other examples of simulation with noise, and other explanations) when 
time < 10, and then continue the simulation without noise up to a time of 500. The stability 
of the (3,2) stationary solution for Neumann boundary conditions is illustrated by the fact 
that the shape returns quickly to this solution (observe the fact that the final amplitude is 
exactly the same as the initial one). 

Of the different stable stationary solutions just described, the largest basin of attraction 
(with initial conditions close to a fiat fiame with some random perturbations) corresponds 
to the most symmetric solution (i.e. the (3,2) solution for 5 poles) and the monocoalescent 
solution ((5,0) in the previous case). It even seems, if one compares both types of solutions, 
that the most symmetric solution has a larger basin of attraction for low values of l/u (in 
the case of five poles for instance), and the monocoalescent one a larger basin for large l/u. 
However, this result could be limited to this type of initial conditions. Actually, in Section 
IIV[ it will be shown that in the presence of a moderate white noise added to the Sivashinsky 
equation, the solution is much more often close to the most symmetric bicoalescent solution 
with the optimal number of poles than close to the corresponding monocoalescent solution. 
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Before closing this section, let us note the analogy of the bicoalescent solutions found 
here with cellular solutions observed experimentally in directional solidification J2|. These 
solutions, called doublets, look almost the same as the bicoalescent solutions of this section. 
They are also stable for some range of parameters. However, a major difference is that 
there is no instability at large scale in directional solidification, and that as a result, the 
structure with one small cusp, one large cusp can be repeated a number of times in the 
overall doublet cellular structure. But in both cases, flames (bicoalescent solutions) and 
directional solidification (doublets), these type of stationary solutions are related to the 
well-known phenomenon of tip-splitting of curved fronts |lfl| . 

III. WEB OF STATIONARY SOLUTIONS 

As most of the solutions of the previous section were not found by Guidi and Marchetti, 
only some trivial, cellular solutions obtained by folding, such as the (2,2) solution, we in- 
vestigate in this section higher values of 1/z/ than those used in their paper (S|. As in this 
paper, we plot the stationary solutions in a diagram giving the amplitude (maximum minus 
minimum value of the solution) versus 1/z/. 

A light version of this diagram, with only the most important solutions, particularly the 
bicoalescent solutions of the previous section, is shown in Figure IH The complete version 
of this diagram, with all the solutions obtained by the author, will be shown in Figure El 
We have found it necessary to use two figures, because the different solutions are so close 
in Figure that it is difficult at first sight to recognize a particular bicoalescent solution 
in this figure. We hope that a comparison between Figures HI and El can help the reader 
understand how the bicoalescent solutions of the previous section are interconnected to the 
rest of the stationary solutions, particularly the cellular ones. But the author knows, it is 
not an easy task for the reader, so for the moment, we only start with the simplified version 
of the diagram. To be more precise, we plot in FigurelHthe basic solutions, i.e. the solutions 
with n poles whose branch exists in the interval [2n — 1, 2?2 + 1] of the parameter 1/z/. In 
this interval, these type of solutions have thus the optimal number of poles, a necessary 
condition for the solution to be stable, as explained in Section ITTl 

In dashed lines in Figure IH can be seen the monocoalescent solutions (n, 0) which are 
created at 1/z/ = 2n — l and are stable in the periodic case until the next solution {n + 1, 0) 
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is created. From these solutions, by a process we call here folding and which is defined 
in the previous section, the solutions (1,1) (1,1,1) ... (dotted lines) are created, as well as 
the bicoalescent (1,1) (2,2) (3,3) ... The non trivial bicoalescent solutions of Section ITTl are 
created starting from these symmetric bicoalescent solutions. The hierarchy (2,1) (3,1) (4,1) 
(5,1)... is created starting from the (1,1) solution obtained by folding. The hierarchy (3,2) 
(4,2) ... emerges from the (2,2) solution. Finally, In FigurelHthe solution (4,3) (first element 
of the hierarchy (5,3) (6,3) ...) is created from the (3,3) solution, which means that one pole 
comes from infinity at a given value of 1/z/ to create the solution. 

All the solutions of the previous hierarchies are plotted as solid black lines in Figure SI 
With the exception of the folded symmetric solutions, all the other bicoalescent solutions 
of this figure are stable when they are created, until a new solution with one more pole 
appears. This behavior is exactly similar to the monocoalescent solutions, the intervals of 
stability are also the same. 

In solid gray lines in Figure IH are plotted however another hierarchy of solutions. This 
hierarchy contains solutions of the type (2,1,1) (3,1,1) (4,1,1) apparently created exactly 
on the same intervals as before. Of course this hierarchy only leads to unstable solutions, 
in the periodic as well as the Neumann case. It seems reasonable to suggest that as 1/z/ 
increases, an infinite number of hierarchies will be created, each starting from a suitable 
folded solution. The author actually suggests the following conjecture: for each value Xjv of 
the control parameter with optimal number of poles n(z/), all the multicoalescent stationary 
solutions having the optimal number of poles, labelled (ni, ...,np) for any -p in the interval 
1 < p < n(z/), with Y3i=\f^i = ^{^)j do exist. 

Furthermore, as the amplitude of the solutions in these hierarchies increases with l/i/, it is 
extremely likely that solutions of the (n, 1) hierarchy for instance, will soon become extremely 
close to the corresponding monocoalescent solution (n + 1, 0). And in the Neumann case, all 
the bicoalescent hierarchies lead to stable stationary solutions. A study of the time evolution 
of solutions of the Sivashinsky equation will be reported in Section IIVI 

The previous argument suggests that there are many stationary solutions of the Sivashin- 
sky equation. However, as seen in Figure Figure El was a very simplified version of the 
diagram, with only the most important stationary solutions, which were called basic solu- 
tions (see the explanation above), and form a sort of skeleton of the entire structure of the 
solutions. We have called this structure web of stationary solutions for obvious reasons, all 
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the solutions are interconnected, even the number of jumps necessary to go from one solu- 
tion to one another can probably be defined, reminiscent of the hops from router to router 
on the internet. It is to be noted that the other well-known Sivashinsky-type equation, the 
Kuramoto-Sivashinsky equation, also admits a huge number of stationary solutions [uj. The 
author does not even claim to have obtained in Figure El something comprehensive in the 
parameter space studied. The reader is again warned that it is easier to look at both figures 
m and ini at the same time, to locate first the basic solutions that a particular interpolating 
solution connects. 

The new solutions compared to Figure IH are of the interpolating type discussed by Guidi 
and Marchetti. We define here these interpolating solutions (as opposed to basic solutions) 
as solutions whose branch does not exist in the interval [2n — l,2n + 1] of the parameter 
1/z/. Thus these solutions do not have the optimal number of poles and cannot be stable 
(starting from such a solution, a pole would come from infinity or disappear at infinity and 
a solution with the optimal number of poles would be created). But in Figure El it can 
be seen that these interpolating solutions typically connect different basic solutions of the 
previous bifurcation diagram (Figure SJ. 

For instance, if one starts from the cellular solutions (1,1,1,...), there exists interpolating 
solutions starting from this solution and leading to all cellular solutions and the mono- 
coalescent solutions above. It must be noted that the precise values of 1/z/, where these 
interpolating branches appear from the cellular solutions, were calculated analytically in 
|l2l |. In the simple case of the (1,1,1) solution already studied by Guidi and Marchetti, it is 
possible to move the poles vertically in the complex plane in two different ways in order to 
have an initial guess of the position of the poles on the interpolating branches (the Newton 
iteration leading to the true values of the positions of the poles). Each interpolating solution 
emanating from a cellular solution can be labelled by the way the poles move along the in- 
terpolating branch compared to the cellular solution. This type of pole movement along the 
interpolating branch (at the beginning, where the branch is created) corresponds exactly to 
the way the poles of the cellular solutions must be moved in order to obtain an initial guess 
that will converge. So we have the (+,-,+) solution: two poles are moved upward in the 
complex plane (i.e. their imaginary part increases, while the real part is kept constant), one 
downward compared to the (1,1,1) solution. This (+,-,+) solution will interpolate, starting 
from the three cells solution, all the monocoalescent solutions (1,0) (2,0) and (3,0) (this part 
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of the diagram will be described in more details later). We have also the (-,+,-) solution, 
which, as seen in the figure, interpolates the (1,1) solution (one pole going at infinity at this 
point). 

If we look at a much more complicated case, the five poles (1,1,1,1,1) solution, it seems 
that in order to get the interpolating solutions, we have to consider at least three levels of 
vertical movement of the imaginary part of the poles, and for instance one interpolating so- 
lution has been constructed by moving the third pole upward, the first and fifth downward, 
the second and fourth somewhere in between. Unfortunately, as shown in the case of the 
interpolating solutions emanating from the six poles cellular (1,1,1,1,1,1) solution, the au- 
thor's capacities have been exceeded and neither the solution interpolating (1,0) (2,0) (3,0) 
(4,0) (5,0) (6,0), neither the one interpolating (1,1,1,1) have been found. Actually, although 
it is more or less obvious that these solutions exist, the present author has been unable to 
generate initial pole locations converging to these solutions (which probably means that the 
author has not understood what type of perturbation of the cellular solution leads to these 
two branches). 

If the way the monocoalescent solutions are interpolated starting from the cellular so- 
lutions is now considered, we prefer to start now from the monocoalescent solution, for 
instance the (6,0) solution, and decrease Xjv. In Figure HJ the monocoalescent solutions 
were appearing suddenly apparently from nothing, for some value of the control parameter. 
On the contrary, in FigureEl precursors of the monocoalescent solution exist. So if the (6,0) 
solution appears at 1/z/ = 11, what do we have exactly before ? 

Actually, between Xjv = 11 and l/u = 10, the precursor of (6,0) is a bicoalescent (5,1) 
solution, with fives poles at zero, one at vr, however, the last one is very far from the real axis, 
and does not lead to a cusp in the solution. This type of bicoalescent solution, apart from the 
folded solutions like (2,2), were the only ones obtained in Guidi and Marchetti (they have 
obtained actually the (3,1) solution interpolating (4,0) and the (2,1) interpolating (3,0)). 
They are unstable even for Neumann boundary conditions, because they do not have the 
optimal number of poles corresponding to the control parameter (the optimal number was 
defined in Section ITT|l . 

Between 1/z/ = 10 and 1/z/ = 9, the solution is no more bicoalescent, but is instead a 
(4,1,1) solution. Then on [8,9] we have a (3,1,1,1) solution, on [7,8] a (2,1,1,1,1) solution, 
and as said before, we have not obtained the precursor close to the six poles cellular solution. 
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It is also possible to explain the previous claim that the precursors of (6,0) interpolate all 
the monocoalescent solutions with a number of poles less than 6. At Xjv = 11, one of the 
six poles at zero goes to infinity, and reappears at vr to give a (5,1) solution. At Xjv = 10 
, the pole at vr and one of the poles at zero go to infinity, and reappear later to give a 
(4,1,1) solution, and so on. The fascinating point is that although all the precursors appear 
different, the curve of the amplitude of all the precursors and of the final monocoalescent 
solution versus the control parameter looks perfecly smooth. This, as well as the overall 
structure of Figure suggests that symmetries less obvious than those leading to the folded 
solutions could be at work in the Sivashinsky equation. 

Let us look now at the shape of all these precursors in physical space. We consider 
as before the case z/ = 10.5 (optimal number of poles : 5) . We show in Figure El different 
curved fiame solutions. The one with the higher amplitude is the stable monocoalescent (5,0) 
solution. Then we have, with smaller amplitude, a six poles (5,1) solution interpolating (6,0). 
Then we have the four poles (4,0) solution, a seven poles (4,1,1,1) solution interpolating (7,0), 
the three poles (3,0) solution, and an eight poles (3,1,1,1,1,1) interpolating (8,0). We have 
stopped there, as the next solutions in this list have an amplitude very different from the 
original (5,0). The interesting point is that in Figure all these solutions, which have a 
very different number of poles, look relatively similar, like subsided versions the original 
monocoalescent solution, the first ones being very close to (5,0) (and will be even closer 
with increasing l/i'). It seems that this is the way the Sivashinsky equation is recovering a 
continuum of curved fiame solutions in the limit l/u -^ oo, something like the continuum of 
Ivantsov parabola of the related solidification problem \\S^. From the simulations of Section 
IIV[ it is not obvious at all that these subsided unstable stationary solutions close to the 
monocoalescent play any particular role in the dynamics, except perhaps by providing ways 
to escape the stable monocoalescent solution during the transient phase. We have shown in 
the successive Figures [ll [HI and IHl the solutions and their poles for the non trivial cases (5,1), 
(4,1,1,1) and (3,1,1,1,1,1) respectively. Once again, it highlights the fact that the presence 
of poles is not equivalent to the presence of cusps, sufficiently far from the real axis, and 
with other poles much closer, some poles only lead to solutions with a weaker amplitude 
(see Appendix EI)- 

Now, if we take another look at the stable bicoalescent solutions of Section Ull the same 
phenomenon as for monocoalescent solutions has to be observed: the bicoalescent solutions 
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do not appear from nothing at a precise value of the parameters, they have precursors, as 
seen in Figured For instance, we have produced precursors of the {n,l) hierarchy, which 
also look like subsided versions of the corresponding stable bicoalescent solutions, and will 
also be closer to the original solution as 1/z/ increases. 

Overall, the bifurcation diagram going from cellular to curved flame fronts with all the 
interpolating solutions of Figure El has a structure totally unexpected. In the Sivashinsky 
equation case, most of the cellular solutions are unstable. However, the addition of a suf- 
ficient amount of gravity (fiames propagating downward) to the Sivashinsky equation is 
known to stabilize these solutions and to create a complex transition from cellular to curved 
.0.S when ,av.. is va,« Q Q- K .e.ai. .o be seen . the s.uctu.e of .Ms t.nsit.n 



has any relation with Figure El which is likely, as a stab 



e stationary solution close to the 
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bicoalescent solutions of the present paper was found in |lj|. But searching for stationary 
solution with gravity is much more difficult than with the Sivashinsky equation, as no pole 
decomposition exists. The author takes this opportunity to say that the instabilities of 
curved fiames observed with a very small gravity (and with periodic boundary conditions) 
in (i3| would probably disappear with Neumann boundary conditions, as the most violent 
instabilities of this paper are created by antisymmetric modes. 

IV. EVOLUTION WITH NOISE 

In Figure ^1 we start by showing a typical time evolution with periodic boundary con- 
ditions, and a white noise added to the right hand side of the Sivashinsky equation. This 
white noise is gaussian, with zero mean value and deviation one, and we multiply it by an 
amplitude a . a = 0.01 and l/u = 11.5 (optimal number of poles: six) in the simulations 
presented in this section, with periodic and Neumann boundary conditions. In Figure HH is 
plotted, for periodic boundary conditions, the amplitude of the front versus time, the initial 
condition being a five poles solution which is not stationary for this value of the control 
parameter, and leads to the initial transient. 

After this transient, the solution oscillates violently between low and high values of the 
amplitude. The peak values correspond to curved front solutions, with the poles being 
apparently almost monocoalescent, but with an amplitude much higher than the monocoa- 
lescent (6,0) stationary solution. The values of the amplitudes for the stationary solutions 



12 



(6,0) (5,1) (4,2) (3,3) are all indicated in the figure by gray lines, so that the reader can 
compare. The low values correspond to shapes with a new cusp formed in the fiat part of the 
front. For the lowest values of the amplitude, this new cusp leads almost to a bicoalescent 
solution, but with again an amplitude which seems higher than the (5,1) or (4,2) stationary 
solution. The solution never comes close to the (3,3) solution, which on [0, 2tt] is a two cells 
solution. Furthermore, other low values of the amplitude correspond to a cusp that develops 
without being exactly centered. Anyway, the dynamics is dominated in the periodic case by 
antisymmetric perturbations. Even if the new cusp formed by the perturbation is correctly 
centered when it forms, it will ultimately move on one side, and will be swallowed by the 
main cusp. This of course modifies the position of the main cusp, and leads to the very 
high peak amplitudes observed. This antisymmetric dynamics is forbidden for Neumann 
boundary conditions, so let us see now what happens in this case. 

The situation is shown in Figure ^2 for the same control parameter and noise amplitude 
as in the periodic case. Before discussing this figure in detail, the overall impression is that 
the signal obtained is much less turbulent. The different stationary solutions for this value 
of the control parameter are also indicated by gray lines. 

The first point to note is that in this figure, except in the initial transient, the front is 
never monocoalescent. Even for the peak values obtained, where the amplitudes obtained 
sometimes seem close to the (6,0) amplitude, we stress that all the solutions obtained at the 
peak value are bicoalescent and not monocoalescent. On the contrary, the solution seems 
often close to the different bicoalescent (3,3) (4,2) and (5,1) solutions. We show in Figure 
[121 a comparison between the solution at time 410.555 in Figure HU (dashed dotted line), 
where the amplitude has a local minimum very close to the amplitude of the (4,2) solution, 
and the shape of the (4,2) solution for l/u = 11.5 (solid line). The agreement between both 
solutions is excellent in this case. For very small values of the noise amplitude (not shown 
here) the solution (with Neumann boundary conditions) actually oscillates around the (3,3) 
solution, without making jumps to any of the other stable stationary solutions. As the noise 
used here is gaussian, it is not impossible however, that jumps could occur as extremely rare 
events (for very small noise amplitudes), and could be observed in very long simulations. 

The value of the noise taken here a = 0.01, although moderate, is already sufficient 
to produce jumps in the amplitude, often actually jumps between the bicoalescent steady 
solutions. The very low values of the amplitude in Figure E] correspond to shapes with three 
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cusps in [0,7r], one on each boundary, and one in the middle. For this value of the control 
parameter, the middle cusp is always smaller than the cusps on the side. It is the author's 
opinion that the lowest values of the amplitude correspond to a shape close to an unstable 
stationary solution, which has not been found in Figure As the solution does not need to 
be symmetric on [0, vr], the mechanism for the disappearance of the middle cusp is relatively 
similar to the same one on [0, 27r] in the periodic case, the middle cusp moves on one side 
and is swallowed by one of the two main cusps. The difference here with the periodic case 
is that the main cusp does not move after having swallowed the small cusp and stays on the 
boundary. 

After the low values of the amplitude comes a transient, where the amplitude very quickly 
grows towards a peak value, which is a very unstationary bicoalescent solution. Depending 
on the noise, the shape will then often come back close to a stationary bicoalescent solution. 
Finally, it seems that higher noise amplitudes or larger Xjv (the type of signal obtained is 
very sensitive to this last value) lead to more turbulent curves of amplitude versus time with 
more jumps and more time spent in the unstable low amplitudes solutions and the very 
unstationary peaks. 

To conclude this section, let us compare the behavior with Neumann and periodic bound- 
ary conditions. For small z/, the stable stationary solutions are very sensitive to external 
noise in both cases. As is well-known in the periodic case (and in this respect, the situation 
is very similar with Neumann boundary conditions) , small perturbations are continuously 
created on the front. But the difference lies in the symmetries. In the periodic case, the sta- 
ble stationary solutions are the monocoalescent solution with the optimal number of poles, 
and the continuum of its lateral translations, all neutrally stable because of this symmetry. 
The noise keeps disturbing the monocoalescent solution, but another solution of the contin- 
uum of monocoalescent solutions (with the optimal number of poles) is also continuously 
recreated. With Neumann boundary conditions, the stable solutions are now the bicoales- 
cent solutions with the optimal number of poles. The perturbations created by the noise 
now serve to explore the different stable bicoalescent solutions, causing jumps between two 
different bicoalescent solutions. But with Neumann boundary conditions, all stable solutions 
are not created equal, some are easier to destabilize than the others. As seen previously 
for instance, the monocoalescent solution is more sensitive to noise. As a result, during the 
time evolution, the front will almost never be close to the monocoalescent solution for small 
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V (which is just the opposite of the behavior with periodic boundary conditions). 

V. CONCLUSION 

To summarize this paper, new bicoalescent solutions of the Sivashinsky equation, stable 
in the Neumann case, have been obtained. They have found their location in the incredible 
structure of the web of stationary solutions. Simulations for moderate noise show that the 
evolution is controlled by jumps between stationary solutions. The author would like to 
insist here on the most important point of this paper: evolution with periodic (controlled 
by antisymmetric perturbations) and Neumann boundary conditions is very different. The 
Neumann boundary conditions are more realistic, although in the presence of heat losses, 
the flame is no more perpendicular to the wall (and is of course three dimensional). Finally, 
it is likely that new analytical studies of the Sivashinsky equation should be possible: even 
if the equation is now almost thirty years old, many things remain to be explained. 

Appendix A: HOW FAR MUST A POLE BE LOCATED FROM THE REAL AXIS 
TO CREATE A NEW CUSP? 

In this appendix, we will try to explain in a very simplified way that adding a new pole 
to a monocoalescent solution does not necessarily create a new cusp if the isolated pole is 
located too far from the real axis. Let us consider the following idealized situation: we have 
a monocoalescent stationary solution with poles located at 0. A new pole at vr is added to 
this solution, without moving any of the other poles coalesced at 0. The front with the new 
pole is no more stationary, but in this appendix, we try to answer the following question: at 
which distance of the new pole to the real axis is a new cusp created ? We call this distance 
yc and its value will be measured numerically for different values of 1/z/, with an optimal 
number of poles. In real situations the presence of the pole at vr modifies the position of the 
poles at 0, particularly the poles located far from the real axis. We neglect this effect as we 
just want to have a reasonable order of magnitude of the value of j/c leading to a new cusp. 

It turns out that the value of yc can be computed analytically in the continuous ap- 
proximation introduced by Thual Frisch and Henon [3|. Instead of summing on every pole 
located at 0, this discrete sum is replaced by an integral, with a pole density p{%j) {y being 
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the vertical coordinate in the complex plane) given by (see [3| for a derivation): 



Q 



n^u \ 4 



Piy) = ^- In ( coth -|- 



The value of the slope of the front ^^^corresponding to the coalesced poles at (in the 
continuous approximation) and to the isolated pole at vr is given by: 

(f)^ (x) = -uP j p (y) cot (— — j dy-u cot ( ^ j - u cot ( ^ j 

where P denotes the principal value of an integral going from — cxd to +cxd, the conjugated 
isolated poles being located at vr ± iyi. As a criterion for the appearance of a new cusp, we 
choose the natural condition (p^x {x = n) < 0. The value of (p^x at this point is created by 
the competition between the coalesced poles at 0, which tend to prevent the creation of the 
new cusp, and the isolated pole (and its complex conjugate) which has the opposite effect. 
With the previous forms of the slope and the pole density, we obtain: 

0..(x = n) = P f^ln (coth IJL\\ J dy-—/ 



27r2 \^ \4: J J cosh^y/2) sinh^ (yi/2) 

Integrating by parts, the antiderivative of the function under the integral sign is 

— In coth — ] I tanh (u) + 2 arctan (exp (u)) j 

with u = y/2 , leading finally to 

4>xx {X = 7l) = - - — -;^ 



TT sinh (1/1/2) 

In this formula, the term I/tt comes from the poles at 0, the other term from the isolated 
pole at n. As said before, these two terms have different signs. The condition (p^x {x = n) = 
finally leads to the value of yi = yc corresponding to the appearance of a cusp, which is: 



yc = 2 arcsinh f ^/txuj 



The cusp only appears if yi < yc- We now compare in Figure UHl this formula to the 
values oi yc measured numerically for l/u =10 (5 poles at 0, one at vr) , 20 (10+1 poles), 
40, 60, 80, 100 (50+1 poles), each time with the optimal number of poles coalesced at 
and one extra pole at vr. The solid curve is the previous formula obtained in the continuous 
approximation, the circles are the values measured numerically. It can be seen that the 
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agreement is good. It is even more surprising if we think that for 1/z/ = 10 we have only 
five poles at and the continuous approximation for the second order derivative at tt works 
correctly, the numerical point is just slightly below the theoretical curve. Of course, this 
result is obtained in the framework of an illustrative model where all the positions of the 
poles are kept fixed, but it serves to justify the fact that in the presence of other poles, a 
new pole at a different x coordinate needs to be sufficiently close to the real axis to create 
a new cusp. 
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Figure 1: Flame shapes (x,0(a;)) with x € [0, vr] of the (from left to right) (5,0) (4,1) and (3,2) 
stationary solutions for 1/u = 10.5. All scales are the same in the x and y direction 



Figure 2: Lower part of the figure (below the horizontal segment) : flame shape (x,0(x)) with 

X € [0, vr] of the (3,2) stationary solution for 1/u = 10.5. Upper part of the figure (above the 

horizontal segment) : corresponding pole locations in the complex plane (the segment is the real 

axis in the complex plane between and 27r, the poles are indicated by circles). All scales are the 

same in the x and y direction, both for the flame shape and for the poles, 
o 



o 

o 



o 
o 
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Figure 3: Amplitude versus time with Neumann boundary conditions for Xjv = 10.5. The initial 
condition is the (3,2) stationary solution. A gaussian white noise (amplitude a = 0.01) is imposed 
on this solution when time < 10, and is then suddenly stopped. The solution goes back to the (3,2) 
solution for large times. 
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Figure 4: Stationary solutions: amplitude A(/) versus \/v (light version with the monocoalescent 
solutions (n,0), the cellular solutions (1,1,1,...) , and the stable bicoalescent solutions (p,q)) 
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Figure 5: Stationary solutions: amplitude A0 versus Ijv (complete version of the solutions obtained 
by the author, including the interpolating solutions). 
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Figure 6: Different curved front solutions (<^(a:)) with x G [0,7r] for \/v = 10.5. A constant has been 
added to each solution in order to have the same spatial mean value for all the solutions presented 
in this figure. 
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Figure 7: Interpolating (5,1) solution for 1/z/ = 10.5: lower part of the figure, flame shape, upper 

part of the figure: pole locations (see Fig. |2lfor a more complete description of this kind of figure) 
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Figure 8: Interpolating (4,1,1,1) solution for 1/z/ = 10.5; lower part of the figure, flame shape, 
upper part of the figure: pole locations (see Fig. [21 for a more complete description of this kind of 

figure) 
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Figure 9: Interpolating (3,1,1,1,1,1) solution for Xjv = 10.5: lower part of the figure, flame shape, 
upper part of the figure: pole locations (see Fig. [21 for a more complete description of this kind of 

figure) 
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Figure 10: Amplitude versus time with periodic boundary conditions for 1/u = 11.5 and a = 0.01 
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Figure 11: Amplitude versus time with Neumann boundary conditions for l/u = 11.5 and a = 0.01 
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Figure 12: Comparison of the solution at time 410.555 (dashed dotted line) in Figure ITTlf Neumann 
boundary conditions l/u = 11.5 and a = 0.01) with the stationary (4,2) solution for 1/v = 11.5 
(solid line) 
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Figure 13: Comparison of the values of yc (maximum distance to the real axis of a pole at vr to 
create a new cusp) obtained numerically with a theoretical value obtained by using the continuous 
approximation of Thual Prisch and Henon 
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